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l/^ Context. With current wide-field near-infrared (NIR) instruments the scattered light in the near-infrared can be mapped over large areas. Below 
A\j ~ 10 mag the surface brightness is directly proportional to the column density, and at slightly higher column densities the saturation of the 
intensity values can be corrected using the ratios of the intensity in different NIR bands. Therefore, NIR scattered light provides a promising 
J new method for the mapping of quiescent interstellar clouds. 

• Aims. We develop a method to convert the observed near-infrared surface brightness into estimates of the column density. We study and quantify 
the effect that different error sources could have on the accuracy of such estimates. We also propose to reduce systematic errors by combining 
surface brightness data with extinction measurements derived from the near-infrared colour excess of background stars. 
O I Methods. Our study is based on a set of three-dimensional magnetohydrodynamic turbulence simulations. Maps of near-infrared scattered light 
are obtained with radiative transfer calculations, and the maps are converted back into column density estimates using the proposed method. 
The results are compared with the true column densities. Extinction measurements are simulated using the same turbulence simulations, and 
^ are used as a complementary column density tracer. 

Results. We find that NIR intensities can be converted into a reliable estimate of the column density in regions with Ay up to almost 20 mag. We 
• • show that the errors can be further reduced with detailed radiative transfer modelling and especially by using the lower resolution information 
_ available through the colour excess data. 

Conclusions. We urge the observers to try this new method out in practice. 

)^ 

^ , Key words. ISM: Structure - ISM: Clouds - Infrared: ISM - dust, extinction - Scattering - Techniques: photometric 



1. Introduction 

The large scale distribution of the interstellar matter is affected 
by Galactic rotation, gravitational instabilities, and supernova 
explosions. Within the clouds the structure results largely from 
supersonic turbulence that is constantly fed by the large scaled 
phenomena and by stellar winds and outflows from newborn 
stars. At smaller scales the self-gravity becomes the dominant 
factor, keeping dense cloud cores together and eventually lead- 
ing to core collapse and the formation of new stars. The mag- 
netic fields are important at all scales. The flow of the matter is 
affected by a number of processes, each having its own imprint 
on the observed density distributions. Investigation of the cloud 
structure is important since it provides clues on the relative im- 
portance of turbulence, self-gravity and magnetic fields at dif- 
ferent scales. It also provides information on different stages of 
the star-formation process and sets constraints on theories of 
star formation. 

Information about the 3D spatial structure of interstellar 
clouds is obtained indirectly, through analysis of the observed 



radiation which represents an integral along the whole line of 
sight. Extended maps of interstellar clouds can be obtained 
with several methods: (1) the integrated intensity of emis- 
sion lines of various molecular or atomic tracers, especially 
CO and HI, (2) the thermal emission of dust grains at far-lR 
and sub-mm wavelengths, (3) star counts at optical or near- 
infrared wavelengths, (3) the near-infrared reddening of the 
light from background stars, (4) mid-infrared absorption to- 
ward dark clouds (Egan et al. IT^^ Hennebelle et al. UMT\ . 
and (5) absorption toward bright x-ray background. Each of 
these methods has limitations, but gives a complementary view 
of the cloud structure. 

A given molecular line is useful only in a limited density 
range, above the critical density of the transition, and below 
the optical depth where the line saturates. Molecular abun- 
dances depend on complicated, time-dependent chemical net- 
works and, furthermore, depletion onto dust grains may cause 
additional effects. In the study of molecular lines, one should 
be able to estimate the excitation of the molecules and possible 
radiative transfer effects. Since the line-of-sight structure of the 
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cloud is unknown, this cannot be done accurately. Abundances 
and excitation vary both on the plane of the sky and along each 
line of sight. It is well known that different molecules may peak 
at entirely different locations and, therefore, the estimated col- 
umn density map depends critically on the selected tracer Low 
extinctions around Ay ~1 pose a special problem because of 
the transition from atomic to molecular gas. Therefore, column 
density must be estimated by combining information from sev- 
eral tracers with largely unknown abundances. The resolution 
obtained by single dish observations is usually some tens of 
arc seconds. Arc second resolution can be reached only with 
interferometric observations of small areas. 

Thermal dust emission at far-IR and sub-mm wavelengths 
could provide a more straightforward picture of column den- 
sity - if dust temperatures can be estimated and the dust to 
gas ratio remains constant. However, both dust temperature and 
FIR/sub-mm optical properties of dust grains may have sig- 
nificant spatial variations (Cambresy et al. '2001' del Burgo et 
al. 12003 Dupac et al. 2003 Kramer et al. 2003 Stepni k et al. 
IMHI Lehtinen et al. (2004 2006 Ridderstad et al. l20n6t . These 
may be due to grain growth by coagulation and ice mantle de- 
position (Ossenkopf & Henning 1994 Krugel & Siebenmorgen 
11994 1 or physical changes in the grain material (e.g., Mennella 
et al. il998i Boudet et al. 12005V Quiescent clouds have rather 
low intensity in the FlR/sub-mm and, because of the lim- 
ited sensitivity, observations have concentrated on regions with 
Ay ~ 10" or above. For single dish observations the spatial 
resolution is typically worse than ~10" in the sub-mm. FIR 
observations must be carried out with space borne instruments 
and, in spite of the shorter wavelengths, the resolution is not 
better 

Optical star counts are used to map the extinction mainly 
at low column densities (Wolf [T923t but can give reliable es- 
timates up to Ay ~5 mag, depending on the available obser- 
vations. High resolution is, however, hard to obtain since each 
resolution element must contain a large number of stars. In the 
optical region the stellar density drops rapidly above Ay ~1, 
but in the NIR the star counting method remains useful beyond 
Ay ~20, provided that deep K-band observations are available. 
However, in NIR the colour excess method yields better spatial 
resolution (e.g., Alves, Lada & Lada 2001 ; 2002). The method 
is statistical in the sense that it relies on average properties 
of the background stars. Unlike star counting, the colour ex- 
cess method provides an independent extinction estimate for 
each star, that is, for a number of very narrow beams through 
the cloud. The errors of these individual extinction estimates 
are dominated by the uncertainty on the (generally unknown) 
spectral type of the background star An actual map of extinc- 
tion is obtained by spatial averaging, and the reliability can 
be improved by combining results from more than two NIR 
bands (Lombardi & Alves 2001) and adaptive spatial resolu- 
tion 1 120021 . With dedicated observations one can reach a res- 
olution of ~10". In the case of the commonly used 2MASS 
survey (limiting Kj magnitude ~15) the spatial resolution is, 
depending on the field location, a few arc minutes and the cov- 
ered extinction range Ay ~ 1-15 mag. 

The scattered light in dark clouds was first detected by pho- 
tographic methods at optical wavelengths by Struve and Elvey 



J1936t and Struve J1936 I). They analyzed the light in terms of 
the dust scattering properties. Later on, Mattila J1970h . ll970b ) 
performed photoelectrical surface brightness observations of 
two dark nebulae and determined the dust albedo and scat- 
tering asymmetry parameter at the UBV bands. Haikala et al. 
( 1995^ made the first imaging of scattered light in a discrete 
diffuse/translucent cloud at the far UV wavelengths. 

The first detection of NIR scattered light in a dark nebula il- 
luminated by the normal interstellar radiation field (ISRF) was 
reported by Lehtinen and Mattila ( 1996 1. They also studied the 
relationship between J, H, and K band surface brightness vs. 
dust column density (measured by NIR colour excess method) 
which they found to be linear up to optical depth of ~ 1 in the 
wavelength band in question (their Fig. 8). From their Monte 
Carlo light scattering calculations they predicted that the linear 
relationship will saturate at optical depths of ~ 1.5 to 2 and then 
turn down. They also made the first determination of the albedo 
of Galactic interstellar grains in the JHK bands. The idea of us- 
ing NIR surface brightness as a high-resolution probe of the 
dust density distribution in optically opaque clouds was pre- 
sented by Lehtinen et al. in the ESO Press release 26a/2003 ' in 
connection with their VLT/ISAAC observations of DC303.8- 
14.2. Nakajima et al. J2003t presented JHK surface brightness 
images of the Lupus 3 dark cloud. They discussed their mea- 
surements in terms of scattered light from dust and presented 
diagrams of JHK surface brightness vs. Ay. 

Padoan, Juvela, & Pelkonen l2006l Droposed scattered near- 
infrared light as a direct measure of the cloud column density. 
That study was motivated by the images of the Perseus region 
obtained by Foster & Goodman (2006), which, in accordance 
with Lehtinen & Mattila ( 1996 1 and Nakajima et al. (2003), il- 
lustrated that large scale mapping of scattered intensity has be- 
come possible even for clouds illuminated by normal ISRF. In 
the near-infrared the dust properties are believed to be rather 
constant, resulting in smaller uncertainties in the scattering 
properties than at shorter wavelengths. In the normal extinction 
curve (Cardelli 1989) the optical depth in the K-band is about 
one tenth of the corresponding value in V-band. Therefore, 
in regions with Ay below 10 mag, the K-band intensity re- 
mains almost directly proportional to the dust column density. 
At higher extinctions the surface brightness starts to saturate, 
the effect being stronger at shorter wavelengths. Padoan et al. 
( 2006 ) showed that if the saturation is taken into account, the 
combination of J-, H-, and K-bands can be used to estimate col- 
umn densities for regions with Ay up to ~20 mag. The method 
avoids many shortcomings of the other methods listed above. 
Most importantly, it provides column density maps of interstel- 
lar clouds at a sub-arcsecond resolution. At low A the method 
is limited mainly by the sensitivity of the NIR observations. 

In this paper we study in more detail the properties of col- 
umn density estimates that are based on near-infrared scatter- 
ing. We use the method presented by Padoan et al. (2006) as 
our starting point. We examine and quantify the effects of pos- 
sible error sources, including differences in the cloud struc- 
ture and column density, dust properties, and radiation field. 
We also address possible complications arising from near- 

' http;//www.eso.org/outreach/press-rel/pr-2003/phot-26-03.html 
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infrared dust emission and diffuse background surface bright- 
ness. Furthermore, we propose to improve the reUabiUty of 
the column density estimates with the help of colour excess 
measurements of background stars (the reddening data are a 
byproduct of the observations) and address the possibility of 
more detailed radiative transfer modelling. 

In Sect.lJlwe present the general method of converting NIR 
surface brightness into column density. In Sect. |3] we present 
the turbulence and the radiative transfer calculations used in 
our tests. The radiative transfer calculations provide simulated 
maps of near-infrared surface brightness that can be converted 
back to column density estimates using the proposed method. 
The results are compared with the true column densities of 
our model clouds in Sect. ^ where we will also examine the 
possibility of using colour excess data of the background stars 
to improve the accuracy of the column density estimates. In 
Sect. |5lwe summarize our conclusions on the accuracy of the 
proposed method and discuss some particular sources of er- 
ror Direct radiative transfer modelling of observations is dis- 
cussed in Appendix |X] and further discussion on some of the 
error sources is provided in AppendixlBl 

2. Conversion of near-infrared surface brightness 
to column density 

Padoan, Juvela & Pelkonen ("2006^ presented a method for 
transforming observations of near-infrared surface brightness 
into estimates of column density. The main points are repeated 
below. We assume that the cloud is illuminated by an isotropic 
radiation field, the observed intensity can be attributed to scat- 
tering from dust particles, and observations of a few near- 
infrared bands are available. Because scattering traces only the 
dust column density, a constant dust to gas ratio is assumed 
when results are transformed into total gas column density. The 
accuracy of these assumptions is discussed later in this paper 

In the optically thin case the intensity of the scattered light 
is directly proportional to the column density. When the optical 
depth becomes close to one, the surface brightness starts to sat- 
urate. In the schematic view of Fig ^ we divide the scattering 
process into two phases. In the first phase the external radia- 
tion enters the cloud and is transported onto the selected sight- 
line where it finally scatters toward the observer In the sec- 
ond phase the radiation is propagated along the selected sight- 
line toward the observer During this phase the intensity is de- 
creased by both absorption and scattering. This schematic view 
is useful to understand the reason why our method is based on 
a formalism (e.g. Eq.[0 that corresponds to the radiation trans- 
fer along an individual line of sight, hence accounts only for 
the second phase. The reason is that the second phase is more 
important than the first for the relation between the surface in- 
tensity and the column density, as the first phase cannot gen- 
erate strong radiation field variations within the cloud tightly 
correlated with the line-of-sight column density. 

Because the incoming radiation field is attenuated during 
the first phase in a way that depends on the density structure 
of the cloud, different locations in the cloud could experience 
different radiation field intensities and this could affect the re- 
lation between the observed intensity and the dust column den- 
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Fig. 1. A schematic view of the creation of the observed sur- 
face brightness of scattered light (see text). The figure shows an 
inhomogeneous cloud with arrows indicating the flow of pho- 
tons. In the first phase external radiation is transported onto 
a selected line of sight where it is scattered toward the ob- 
server. Photons reach this line possibly after several scatterings 
and preferably through regions of low density. In the second 
phase radiation propagates out from the cloud along the se- 
lected sightline. 



sity, but this effect is relatively small for two reasons. First, a 
volume element is illuminated from all directions. In an inho- 
mogeneous cloud the radiation propagates more freely through 
low density regions so that inside the cloud the field can be 
much stronger than what would be expected based on the mean 
optical depth. Second, photons can reach a given line of sight 
also after some scatterings. Most scatterings are in the forward 
direction so the radiation can again better reach deeper cloud 
regions. Both effects make the field strength resulting from the 
first phase relatively uniform within the cloud. During the sec- 
ond phase, instead, the attenuation depends on the optical depth 
(column density) of a particular line of sight and every absorp- 
tion and scattering event decreases the observed intensity. 

In the near-infrared the albedo of dust grains is close to one 
half (Lehtinen & Mattila 1996 1, so the scattering is a signifi- 
cant component of the total extinction. Taking into account the 
cloud inhomogeneity, the effect of the attenuation during the 
first phase is perhaps half of the effect during the second phase 
and, therefore, gives a smaller contribution to the relation be- 
tween surface brightness and line-of-sight column density. The 
uncertainty in our method due to this small contribution from 
the first phase is studied with radiative transfer simulations. 

The relation between surface brightness and column den- 
sity follows from the solution of the radiative transfer prob- 
lem. As explained above, the observed intensity is determined 
mainly by radiative transfer effects during the second phase. 
This suggests that a radiative transfer equation written for an 
individual line of sight should provide a good functional form 
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for the solution. In the case of a homogeneous medium this 
results in the equation 

/, ^fld-e-"^), (1) 

which gives the surface brightness of scattered light, ly, as a 
function of column density, A^. In this formula, a and b are con- 
stants related to the dust properties and to the incoming radi- 
ation field. The constant b is clearly similar to an absorption 
coefficient. This analytical formula was already found to fit the 
simulation results of Padoan et al. (2006). We will examine 
these correlations for a larger set of cloud models in Sect. 14. II 
In the limit of low optical depth Eq.[T]gives 

I^abN, (2) 

the intensity is directly proportional to the column density, and 
the product a x b describes the scattering that takes place per 
unit column density. If b is assumed to describe dust properties, 
the constant a would mostly reflect the intensity of the incom- 
ing radiation. However, the interpretation of these constants is 
not this straightforward. In the term exp(-bN), the constant b 
has the role of an extinction coefficient, while in Eq.|2it would 
be a scattering coefficient. In the following we treat a and b 
simply as empirical parameters. We expect that if either the 
dust properties or the radiation field were changed, the values 
of both a and b would change. 

At large optical depths one must correct for the non- 
linearity, i.e., the saturation of the surface brightness. Eq^ 
written for each band separately, defines a parametric curve in 
the (/j. III, /K)-space. Each observed (/j. In, /K)-triplet should 
correspond to a point on this curve, and the parameter value 
could be calculated if the constants a and b were known. We 
emphasize that Eq. simply defines an empirical curve rep- 
resenting the observations. It gives a good description of the 
relation between surface brightness and column density, but is 
not necessarily the optimal function. 

When no independent column density estimates are avail- 
able, surface brightness observations can only be used to deter- 
mine the ratios between the NIR bands. For example, by writ- 
ing Eq.^for the H and K bands, can be eliminated and the 
H band can be expressed as a function of the K band, 

/H = flHX(l-(l--^)^). (3) 

The constants oj, an, and ok and the ratios Zjj/Zjk and bnjb-^ 
can be determined by fitting this curve to the observations. 
However, determination of the individual b constants requires 
additional information. 

We have explained above that the b constants depend 
mainly on the properties of the dust grains. In the NIR, dust 
properties vary only relatively little. If the b constants are sim- 
ilar in different clouds, one can use values already established 
from observations of other similar objects. The a constants de- 
pend mainly on the radiation field that illuminates the cloud. 
The uncertainty of the intensity is often large, ~50%. This af- 
fects mostly the absolute scaling of the column density map. 
The spectral shape of the incoming radiation can also vary de- 
pending, for example, on the presence of nearby OB associ- 
ations. However, if the relative values of the scattering cross 



sections are known, the spectrum can be estimated using ob- 
servations of low extinction sightlines. 

The expected values of the constants a and b can be de- 
termined by radiative transfer modelling, while the observed 
ratios of Eq.|3lcan be used to correct the values of a and the ra- 
tios between the b constants. If a correction is necessary, either 
the radiation field or the NIR dust properties differ from their 
initially assumed values. Cloud properties, the average column 
density and the inhomogeneity, are also expected to have some 
effect. These will be studied in Sect.|3l 

There is still another, more straightforward way to check 
the validity of Eq. ^ When NIR scattering is measured, the 
observations contain colour excess data for a large number of 
background stars. The resolution of the extinction map may be 
one or two orders of magnitude worse than the resolution of the 
surface brightness maps. However, it is sufficient to check the 
parameters of the surface brightness method and even to de- 
termine the values of the individual b constants. It is not clear 
whether the colour excess method is intrinsically more accu- 
rate. However, since the error sources and error properties of 
the two methods are very different, a comparison should be use- 
ful to test for systematic effects and to provide error estimates 
for the derived column density maps. 

Once the constants are determined, Eq. ^ can be used to 
convert surface brightness values into column density. The con- 
version could be done for each band separately, by inverting 
Eq.C] 

N^-Uogil-hla), (4) 
b 

In principle, observations at one wavelength would be suffi- 
cient. However, the use of several bands is more secure, and it 
brings additional information that can be used in the analysis 
(see, e.g.. Appendix. |AJ. In this paper we assume that obser- 
vations exist of J-, H-, and K-bands. Because of the noise in 
the observations and possible model errors the observed inten- 
sities do not fall exactly on the curve given by Eq.^ For each 
map position, we find on this curve a point that minimizes a 
least square distance from the observed (/j. In, /K)-triplet. The 
corresponding column density is then obtained from Eq.|3 

3. Modelling of NIR scattering 

In this section we describe the numerical modelling used 
to study the relation between the NIR surface brightness 
and the column density. The results are based on three- 
dimensional magnetohydrodynamic (MHD) turbulence simu- 
lations and Monte Carlo radiative transfer calculations. 

3.1. MHD turbulence simulations 

The density structure of our model clouds is based on six MHD 
turbulence simulations that provide a good approximation of 
the structure of interstellar clouds. 

The first three model clouds - A, B, and C - have a res- 
olution of 128^ computational cells. In these models the tur- 
bulence is highly supersonic, with sonic rms mach numbers 
Ms ~10. The initial density and magnetic fields are uniform. 
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Fig. 2. Examples of the structure seen in the model clouds. 
The frames a - d show column density maps of models A, C, 
D, and E, respectively. The colour scales are logarithmic and 
independent for each frame. 



and an external random force is applied to drive the turbu- 
lence at a roughly constant rms Mach number In model C self- 
gravity is also included. Model A has approximate equipartition 
of magnetic and kinetic energy, the Alfvenic Mach number be- 
ing Myi ~1. The models B and C have weaker magnetic fields 
and the turbulence is strongly super- Alfvenic with Ma ~10. 
These MHD simulations are the same used in Juvela, Padoan 
& Nordlund ( BOOl ) for estimation of molecular line cooling in 
inhomogeneous clouds, and the reader is referred to that article 
for further details. 

The other three model clouds - D, E, and F - are from 
MHD simulations on grids of 250^ computational cells, but the 
models are downsized to 125^ cells for the radiative transfer 
calculations. The rms sonic Mach numbers are 0.6, 2.5, and 
10.0, for D, E, and F, respectively. The initial density and mag- 
netic field are uniform and the gas is assumed to be isothermal. 
The initial Alfvenic Mach numbers are Ma ~10. The volume- 
averaged magnetic field strength is constant in time because of 
the imposed flux conservation. The magnetic energy is instead 
amplified. The initial value of the ratio of average magnetic 
and dynamic pressures is (Pm)in/(^'d)in = 0.005, so the run 
is initially super- Alfvenic. The value of the same ratio at later 
times is larger, due to the magnetic energy amplification, but 
still significantly lower than unity, (/'m)in/(^d)in = 0.12. The 
turbulence is therefore super-Alfvenic at all times. Turbulence 
is again set up as an initial large-scale random and solenoidal 
velocity field and is maintained with an external large-scale 
random and solenoidal force. Experiments are run for approxi- 
mately 10 dynamical times in order to achieve a statistically re- 
laxed state. The models D-E were used, for example, in Juvela, 
Padoan & Jimenez (I2003> for estimation of [CII] cooling of in- 



homogeneous translucent clouds. Column density maps from 
models A, C, D, and E are shown in Fig.|2l 

Supersonic and super-Alfvenic turbulence of an isothermal 
gas generates a density distribution with a very strong con- 
trast of several orders of magnitude. It has been shown to pro- 
vide a good description of the dynamics of molecular clouds 
and of their highly fragmented nature (e.g., Padoan et al. 1999) 
I2001li2004t . The most important property affecting the radia- 
tive transfer calculations is the degree of inhomogeneity. In this 
respect the models provide a wide range of conditions. Details 
about the numerical method used in the MHD simulations are 
given in Padoan & Nordlund t' 1999V All simulations use peri- 
odic boundary conditions, the effect of which is examined later 
in Sect.lO 

3.2. Radiative transfer calculations 

The flux of scattered radiation is calculated with a Monte Carlo 
program (Juvela & Padoan 2003 , Juvela 2005 ), where sampling 
of scattered radiation is further improved with the 'peel-off' 
method (Yusef et al. 1994 1. During each run, photons are sim- 
ulated at one wavelength, and the scattered intensity, including 
multiple scatterings, is registered toward a selected direction. 
The result is an intensity map where the pixel size corresponds 
to the cell size of the model cloud. The maps size is therefore 
either 125x125 or 128 x 128 pixels. Maps are calculated for the 
three NIR bands J, H, and K (1.25, 1.65, and 2.2 yum), for three 
directions perpendicular to the faces of the cubic model cloud 
(directions X, Y, and Z) and for two diagonal directions (Dl 
and D2). The final maps are averages of several independent 
runs, allowing us to control the level of the numerical errors. 
In the following the uncertainties are characterized with rms 
values, 

rms(x) = ^Yj^j/N, (5) 

where x, form a sample of independent observations of the 
variable x. In the simulations the Monte Carlo noise, i.e., the 
rms-value of A/// resulting from the simulation procedure it- 
self is in all cases below 2%. Because of the moderate optical 
depths and the use of both forced first scattering and peel-off 
methods, the errors are relatively uniform, the relative uncer- 
tainty still being somewhat higher in regions of low column 
density. In simulations of actual observations, additional noise 
is added to these maps in order to simulate the measurement 
errors. At low intensities the added noise component is always 
large compared with the Monte Carlo noise. At high intensities 
the Monte Carlo noise becomes dominant, but is still below 
~2% when the surface brightness exceeds the average intensity 
of the map. 

In the radiative transfer calculations the model clouds can 
be scaled arbitrarily, and the obtained surface brightness is af- 
fected only by the column density. As an example, we scale 
all the model clouds to a size of L = 1.0 pc, and rescale the 
densities to obtain model clouds of different opacity. In the 
case of the 125^ cell models D-F one cell would correspond 
to 0.004 pc (825 AU) and, if that cloud were at the distance of 
415 pc, the pixel size would be one arc second. 
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Fig. 3. Histograms of Ay values in the six model clouds with 
average visual extinction < Ay >= 1.6 mag. For model C the 
distributions are plotted for three orthogonal directions of ob- 
servations. For the other models distributions are shown only 
for one direction. 



Initially we assume a mean density of lO-'cm"^. This re- 
sults in an average column density of 3.1x10^' cm"^ which, 
in the case of the normal Milky Way dust (R = 3.1), corre- 
sponds to 1.6 mag of visual extinction. The actual range of ex- 
tinctions depends on the inhomogeneity of the models. Fig |3l 
shows the Ay distributions for the six model clouds in the case 
of < n >= 10^ cm"^. For example, in the model D practically 
all lines of sight have extinction below 3 magnitudes, while 
in the model C the maximum extinction is over 17™. Later 
we consider also some models with twice as high density and 
< Ay >= 3.2 mag. In the figures we will use visual extinction 
instead of column density. Because the dust properties are con- 
stant within each model, the ratio between Ay and true column 
density is also always constant. 

In the simulations the model cloud is illuminated by an 
isotropic radiation field, with intensities calculated according to 
Mathis et al. ( 1983 ) model of the ISRF. However, since the ob- 
served surface brightness is directly proportional to the inten- 
sity of the incoming radiation, the simulations can be re-scaled 
for any spectrum of the background radiation. Based on DIRBE 
data Lehtinen & Mattila J1996> estimated the near-infrared in- 
tensity of the ISRF to be higher by some 50%. In our analysis 
this higher intensity level would only affect the estimated ob- 
serving times, which would be shorter by a factor of ~2. 




Fig. 4. A map of the difference between the simulated surface 
brightness and the prediction based on a curve fitted according 
to Eq. n The cloud is model C with < Ay >= 1.6. The fig- 
ure shows the shadowing effects that are caused by regions of 
moderate optical depth. 

Dust properties are based on Draine (I2003> and we use the 
data files available on the web ^. In most cases we assume nor- 
mal MiUcy Way dust for which the ratio of total to selective 
extinction is Ry - Ay/E(B -V) - 3.1. For the scattering func- 
tion we use the tabulated scattering phase functions that are 
available on the web. In Sect. I4.2.2l we will examine the effect 
of dust with Ry = 4.0 or Ry = 5.5, and the effect of spatially 
varying dust properties. In those cases the scattering calcula- 
tions are based on the Henyey-Greenstein scattering functions 
(Henyey & Greenstein lI941> and the asymmetry parameters g 
calculated for the different Ry cases. 

4. Results 

In this section we examine how the correlation between cloud 
column density and the observed NIR surface brightness de- 
pends on factors like the cloud structure, the average opti- 
cal depth and the dust properties. The methods of Sect. |2] are 
used to convert surface brightness measurements into column 
density estimates. We study the accuracy of the obtained col- 
umn density maps and examine how the results could be im- 
proved with detailed radiative transfer modelling or by the use 
of colour excess data from background stars. 

4.1. Correlation between scattered light and column 
density 

If NIR scattered light is to be used for a reliable estimation of 
the column density, the relation between these quantities should 
be universal or, when differences occur, those variations should 
be taken into account using the available observations. In the 
following we study how the results are affected by the cloud 
properties. 

^ |http : //www ■ astro . princeton . edu/ draine/dust/| 
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Fig. 5. Relative rms-variation between between the simulated 
//-band surface brightness and the average curve based on 
Eq. n] for models A, C, D, and E, where the mean visual ex- 
tinction is scaled to either A = 1.6 (thin lines) or A - 3.2 (thick 
lines). For each model separate curves are shown correspond- 
ing to observations from three different directions {X, Y, and 
Dl). This plot illustrates the amount of variation in the surface 
brightness of scattered light for a given line-of-sight column 
density. 



Figure |3 shows the difference between simulated surface 
brightness and predictions of Eq. that best fits the points 
(NJviH)). The differences are dominated by systematic fea- 
tures that are connected with the cloud structure (see Fig.lJJi) 
and arise from the attenuation of the radiation field before 
reaching the line of sight in question. This is most noticeable 
in the case of thick filaments. Toward the cloud centre the fil- 
aments appear darker than on the outer side that is subjected 
to a stronger radiation field. When the method of Sect. |2l is 
applied, this effect will probably dominate the pixel-to-pixel 
errors in the derived column density maps. In Fig. |5] the rela- 
tive rms-errors, nns(AIy/Iy), are plotted for selected models as 
functions of the true Ay values. 

A second source of uncertainty is the fact that parameters 
of Eq.n might vary from source to source in an unpredictable 
fashion. This could produce more systematic errors, such as a 
bias in all column density values or in the ratio between low 
and high column densities. We will later study to what extent 
such variations can be corrected using the observations them- 
selves. Here we consider only the amount of variation caused 
by the Ay and the cloud structure. Figure |6l shows for some 
models the fitted curves of Eq.Q] The deviation from the aver- 
age behavior is largest for the most optically thick model C. In 
this case the maximum extinction is high, over 30 magnitudes. 
The shadowing produced by such optically very thick regions is 
visible even at lower Ay as relatively lower surface brightness. 

In Fig.0shows scatter plots of surface brightness versus Ay 
for data combined from all six model clouds and three viewing 
directions (X, Y, and Dl). The two frames correspond to the 
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Fig. 6. Curves of Eq. ^ fitted to data from selected models 
(models as in Fig.|5}- The lower curve is for the K band and the 
upper curve for the J band. Each curve is drawn for the range 
of Ay values found in the corresponding model and direction 
of observations. 




Av Av 



Fig. 7. Scatter plots of surface brightness versus Ay. The plot 
combines data from all six model clouds and three directions 
of observation (X, Y, and Dl). The two frames correspond to 
average visual extinction < Ay >=1.6 and 3.2. The colour scale 
is linear with respect to the density of points. 



two different values of average extinction. The figure shows the 
total expected dispersion between surface brightness and Ay. It 
therefore illustrates the accuracy of a column density estimate 
based on the surface brightness of NIR scattered light. 

In real clouds high densities are limited to the clouds inner 
parts and are, therefore, never subjected to full, unattenuated 
external field. On the other hand, the MHD turbulence sim- 
ulations employ periodic boundary conditions and large den- 
sities may be found near the "surface" of our model clouds. 
Furthermore, when we look at a model perpendicular to its 
sides (direction X, Y, or Z), we have sightlines with high den- 
sities on both the front and back side of the cloud. This means 
that in our simulations the scatter in the relation between sur- 
face brightness and column density is larger than in real clouds. 
We examined the magnitude of this effect for all the models 
with < Ay >= 3.2™. We masked out five-pixel-wide map bor- 
ders and sightlines where the density on the front surface ex- 
ceeded one third of the average density of the cloud. The scatter 
with respect to the analytical fit of Eq.Qldecreased by less than 
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~10%, showing that edge effects will have only a small influ- 
ence on the results. 

4.2. Accuracy of the estimated column density 

In this section we use the method of Sect.|2]to convert simulated 
surface brightness maps into estimates of column density. We 
examine the effects that column density distribution and dust 
properties can have on such estimates. 

There are several possible ways of using Eq. ^ One could 
use either a generic set of parameters a and b, or try to im- 
prove the parameter values using the observations themselves. 
For example, the correct ratio of a-parameters can be deter- 
mined directly from observations at low column densities. In 
practice such a procedure might be necessary, because, for ex- 
ample, the spectrum of the illuminating radiation field might 
not be known beforehand with sufficient accuracy. These pos- 
sibilities are examined further in Sect. l4.3l and in AppendixIXl 
In this section we use a single set of values of the a and b con- 
stants that are obtained as the average of all model runs. As 
shown by Fig.|6l there is little variation between the different 
models, except for the model C with < Ay >=3.2, where the 
presence of optically thick knots has changed the relation even 
for sightlines with Ay ~10. This extreme model was excluded 
from the calculation of the average parameters. 

We will also consider the effect of random noise to simulate 
the effect of observational errors, assuming signal to noise ra- 
tios of 20, 15, and 7 for the J, H, and K bands. We will take the 
ISAAC/VLT instrument as an example. Assuming a 0.6" res- 
olution the total integration time (on-source) is approximately 
50 hours. In the following we assume that one pixel in the sim- 
ulated maps corresponds to 0.3". The S/N ratio per pixel is 
lower by a factor of two, but in the subsequent tests we smooth 
the maps down to the resolution of two pixels. Consequently, in 
the tests the actual noise corresponds to the S/N ratios quoted 
above. The ratios were calculated for the expected average sur- 
face brightness of a cloud with < Ay >= 1.6. In the denser 
regions and in the more optically thick clouds the relative noise 
will be significantly lower Compared with Padoan et al. 2006 
our integration times are longer, mainly because of the lower 
average column density of the models. 

The accuracies derived in this section reflect the intrinsic 
variation between clouds with different density structure and 
different optical depth, and the effect of observational eiTors. 

4.2.1 . Effect of column density distribution 

Fig.Elshows some maps of the error in the column density es- 
timates when no observational noise is added. The input maps 
were convolved with a beam with fwhm equal to twice the pixel 
size so that the estimated Monte Carlo noise is brought down 
close to one per cent. Because the S/N ratio is the same for 
all bands, the accuracy is limited mostly by the saturation of 
the surface brightness and can be improved by giving a larger 
weight to the longer wavelength data. In this case the optimal 
ratio was 1:3:12 for the weighting of J-, H-, and K-band data. 
In the case of actual observations, the different S/N ratios will 



change the optimal ratio. One would also expect the overall 
rms error to be reduced if the weight of the shorter wavelengths 
were decreased according to the extinction. In practice this was 
found to yield only little improvement, and the column density 
predictions were calculated using the constant weighting ratio. 
The errors were largest in the case of model C, where, after the 
smoothing of the maps, the maximum visual extinction was 14 
magnitudes. The maximum error was 23%. The average rela- 
tive rms error was only 1.7%, large part of this number still 
being noise from the simulation procedure. This is to be ex- 
pected, since the majority of sightlines are optically thin and 
the surface brightness is directly proportional to the column 
density. In model C the largest deviations are found in dense 
filaments. The other clouds are more homogeneous, and the 
errors reflect large scale gradients in the strength of the radia- 
tion field: The column density is systematically over-estimated 
close to the cloud borders. 

In Fig. Owe show the average rms-error of the column den- 
sity estimates as a function of Ay. As before, the predictions 
are calculated using the average values of the parameters a and 
b. In models with < Ay >=1.6 the errors remain mostly below 
10% with the exception of model C where, as akeady men- 
tioned, the largest error is ~23%. In the optically thicker clouds, 

< Ay >= 3.2™ , the errors are larger for all sightlines, as the 
shadowing by optically thick regions causes fluctuations in the 
strength of the radiation field. However, below Ay ~ 10™ the 
errors are at most 10%. Above Ay ~ 20™ the average errors 
increase close to 50%. 

Next, we add noise to the surface brightness maps to simu- 
late the observational uncertainties. The S/N ratios are 20, 15, 
and 7 for the J, H, and K bands respectively, when calculated 
for the average surface brightness of the < Ay >=1.6 maps 
and the maps are smoothed to a resolution of two pixels. The 
weighting of the different bands is kept the same as before. The 
results are shown in Fig.^| The main effect of the added noise 
is a decrease in the accuracy at small column densities. The rel- 
ative errors are ~IQ% at Ay = 0.5™, decrease until A ~ 4™, and 
again increase at higher column densities as in Fig.|9l A good 
signal-to-noise ratio should be particularly important at high 
opacities where, due to the saturation of the surface brightness 
values, a small change in surface brightness corresponds to a 
large change in column density. However, in these regions the 
S/N ratios are already very high and, beyond Ay ~ 8™, Fig.|9] 
and Fig. [lO| are practically identical. Figure ^2 shows maps 
of the estimated column density maps for the four models of 

Fig.m 

The column density distribution also depends on the res- 
olution of the observations or, in our case, the resolution of 
the simulations. The discretization of the cloud models could 
affect the results if individual cells were optically thick and 
if a finer discretization resolved significant density variations 
within the original cells. The largest density contrast is found 
in model C. When the cloud has an average visual extinction of 

< Ay >= 3.2™, the 7-band optical depth of the densest cell is 
~1.3 and 99.99% of the cells have optical depth below rj=0.5. 
The maximum optical depth is smaller for the other bands and 
models and, naturally, for the clouds scaled to < Ay >= 1.6. 
Therefore, discretization is not expected to produce errors that 
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Fig. 8. Maps of relative error of column density estimates in the 
case of the models shown in Fig. |2] The average extinction is 
1.6 magnitudes (see Fig.|3, the input maps are convolved with 
a beam with fwhm equal to two pixels, and no observational 
noise is added. The contours are drawn at the levels of +5% 
and +10%. The plots show the total projected area of model 
clouds. 



could affect our results. This was still checked using the model 
F. The original MHD simulation was done on a grid of 256^ 
cells and based on that we constructed three models consist- 
ing of 256"^, 125^, and cells. Radiative transfer simulations 
were repeated for all the three model clouds, and Fig.ll2lshows 
the errors when column densities are estimated at the full res- 
olution of each model. The larger differences at low Ay are 
caused by a difference in the Monte Carlo noise of the simula- 
tions. Otherwise, the accuracy of the column density determi- 
nation is not seriously affected by the numerical resolution. 

4.2.2. Effect of radiation field and dust properties 

In deriving the column densities we have so far assumed that 
model parameters a and b have already been determined from 
other models illuminated with similar external radiation field 
and containing dust with exactly the same scattering properties. 
In reality, both radiation field and dust properties can be subject 
to significant variations. 

If the total level of the incoming radiation is underesti- 
mated, the column density is correspondingly overestimated. 
Below Ay ~ 10" the surface brightness is directly propor- 
tional to the external field so that an error in the assumed in- 
tensity of the radiation field results in similar relative error in 
the column densities. This is a constant factor and does not af- 
fect relative accuracy between map positions. An error in the 
assumed spectrum, i.e., in the intensity ratio between differ- 
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Fig. 9. Rms-error of column density estimates as a function of 
the line-of-sight visual extinction. The plot includes results for 
all six model clouds, for three different directions of observa- 
tion (X, Y, and Dl). The frames correspond to the two values 
of average visual extinction. 



ent bands, can result in different kinds of errors, especially if 
the bands are weighted differently, according to the local S/N- 
ratios. However, observations of optically thin sightlines can 
be used to correct for changes in the spectrum and, partially, 
changes in dust properties. 

We repeated simulations of model C using dust with Ry - 
4.0 and Ry -5.5. The dust parameters were taken from Draine 
et al. ('20031. The J- and H-band intensities were plotted against 
the K-band intensity, and Eq.|3was fitted to the data. The pa- 
rameters flK and Zjk were kept constant using the previous val- 
ues that were derived assuming /?v = 3.1. In the NIR the shape 
of the extinction curve remains rather constant, and the main 
effect results from differences in the absolute values of the 
absorption and scattering cross sections. In the optically thin 
regime the product axbis the slope of the relation between sur- 
face brightness and column density. Since these ratios of a xb 
in different bands can be determined from the data itself, no 
a priori assumption needs to be made about the spectrum of 
the external radiation, and some of the uncertainty caused by 
unknown dust scattering properties is eliminated. However, be- 
cause the value of the product ok x b^ must be assumed, the 
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Fig. 10. Rms-error of column density estimates as a function 
of A. The plot includes results for all six model clouds and 
for three different directions of observation (X, Y, and Dl). 
Observational noise was added to the input maps (see text). 
The frames correspond to the two values of average visual ex- 
tinction. 

column density estimates may be at error by a constant factor 
as discussed above. The situation is worse at high column den- 
sities because we cannot directly determine the column density 
at which the saturation of the surface brightness values starts. 
While the ratios between the b parameters can be recovered 
from the data, the value of an individual b parameter must be 
assumed a priori. If the dust extinction cross sections are dif- 
ferent from the assumed one, the saturation begins at a higher 
or lower column density than expected, and all column den- 
sity estimates of high extinction sightlines are correspondingly 
erroneous. Moreover, if the radiation field is stronger than ex- 
pected, the column densities will be overestimated in the non- 
linear part of Eq. ^ If the surface brightness exceeds the ex- 
pected asymptotic value, represented by the parameters a, no 
solution can be found. However, in such a case the peak sur- 
face brightness can be used to correct the value of ak. 

Figure[Olshows the bias and statistical noise in the cases of 
Ry = 3.1 01 Ry - 5.5. The b^ and parameters correspond 
to a case Ry - 4.0, and the other parameters have been re- 
estimated from the plots of the J- and H-band intensity versus 
the K-band intensity. In frame a the radiation field is the same 
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Fig. 11. Estimated column density maps of the four models 
(A, C, D, and E) of Fig.|3 The column densities are calculated 
based on surface brightness maps where observational noise 
is included (see text). The images show the logarithm of the 
column density, transformed into Ay. The frames a, b, and d 
share a common colour scale. The model D has a much smaller 
range of column densities (see Fig. O, and the frame c has a 
different colour scale. 




3 4 5 

Ay [inag] 

Fig. 12. The errors of the estimated column densities when 
model F is discretized into 64^, 125^, or 256^ cells. The esti- 
mates are based on average parameters derived from all models 
with < Ay >- 1.6™ and < Ay >- 3.2™. The points show the 
average bias in different Ay intervals and the errorbars corre- 
spond to the 1-cr variation within each interval. At low Ay part 
of the differences is caused by the Monte Carlo noise of the 
simulations, which was higher for models of higher resolution. 
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as before. No additional observational noise was included. The 
surface brightness of the scattered light increases with the 
value of the dust. Therefore, for Ry = 5.5 the column densi- 
ties are overestimated and on the non-linear part the systematic 
error increases rapidly as the line-of-sight extinction exceeds 
Ay ~ 10™. The maximum observed surface brightness was 
not used to correct the value of the parameter qk- In the case 
Rv - 3.1 the column densities are underestimated, but errors 
are smaller, being ~2Q% or less for Ay ~ 10™. 

In frame b the radiation field was increased by 30% rel- 
ative to the value assumed in the column density derivation. 
At higher column densities, and especially in the case with 

- 5.5, the column densities would be either seriously over- 
estimated or no value could be obtained from Eq.[2 Therefore, 
the analysis was done using a corrected ak value that was set 
to the maximum observed K-band intensity. Such a correction 
assumes that maps contain at least one optically thick region. 
Because we have assumed that the radiation field is underesti- 
mated, the column densities are overestimated for all sightlines 
and for both dust models. However, in the case of dust with 
Ry - 5.5 the errors are, in fact, much smaller than before. 
Apparently the correction of the parameter ak compensated 
also for the higher than expected dust scattering efficiency. On 
the other hand, for dust with Ry - 3.1 the corrected ak value 
does not improve the estimates. Without clearly saturated re- 
gions the maximum surface brightness is a poor measure of ok- 
In that case the parameter values can be corrected only by other 
methods, for example by comparison with an extinction map 
(see Sect. 14. 3> . In some cases the accuracy can be improved by 
more detailed modelling of the observed source. This possibil- 
ity is discussed further in AppendixlAl 

4.3. Use of extinction measurements 

Observations of NIR surface brightness will automatically give 
photometry for a large number of background stars. The colour 
excesses of the stars can be converted into estimates of the 
clouds column density (Alves, Lada & Lada .2001 , Lombardi 
& Alves l200n Cambresy et al. 2002 1. Each star can be treated 
as a probe of the extinction along individual, very narrow sight- 
lines. Usually the results are spatially smoothed to give a map 
of extinction at lower resolution. Averaging is required because 
of the scatter in the intrinsic colours of the stars. The resolu- 
tion of the resulting column density map is much lower than 
the resolution of the surface brightness maps. This information 
can, however, be used to assist the use of surface brightness 
measurements. The combination of the two methods is inter- 
esting because, while both rely on absorption and scattering 
by dust particles, their actual dependence on the dust proper- 
ties is different. Furthermore, the colour excess method is in- 
dependent of the local radiation field. This makes it possible 
to take into account any radiation anisotropics in the studied 
field. Conversely, differences in the two column density esti- 
mates serve as a sensitive indicator of variations in dust prop- 
erties and the radiation field. 

In the following we use simplified simulations of the colour 
excess method where extinction maps are calculated with the 
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Fig. 13. Effect of a change in dust properties and/or field 
strength on the column density estimates. The plots show the 
bias and scatter in the estimates when dust properties are as- 
sumed to correspond to Ry - 4.0, while they in reality corre- 
spond to either Ry - 3A or Ry = 5.5. In frame b the radiation 
field is also 30% higher than assumed in the column density 
derivation, and the parameter a/j- is re-estimated based on the 
maximum surface brightness. The small circles show the result 
in the case of Ry = 3.1 without re-evaluation of parameter ok- 

NICER method (Lombardi & Alves 1200 II. Stars are placed at 
random locations behind the model cloud and their extinction 
is proportional to the true column density. The magnitude dis- 
tribution follows a linear (mag, logN) relation with a slope of 
0.35. A random error of 0. 12 magnitudes is added to each band 
in order to simulate photometric errors and the variation in the 
intrinsic colours of the stars. The resulting dispersion in the 
colours is similar to or slightly smaller than the scatter found 
in 2MASS data. For simplicity, the same noise was applied to 
all stars. Faint stars are removed when their magnitudes exceed 
some limit rriK in the K-band, niK + 0.85 magnitudes in the 
H-band, or nik + 1.5 magnitudes in the J-band. The iuk limit 
determines the remaining number of stars. 

The comparison of extinction and surface brightness gives 
a way to determine all the parameters of Eq. [2 directly from 
observations. This was tested on model C, using the 128x128 
maps and adding observational noise as in the previous section. 
The upper frame of Fig.^^shows the accuracy of the extinction 
maps that were derived from simulated colour excess observa- 
tions and were averaged with a Gaussian with FWHM equal to 
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7 pixels. The number of background stars over the whole map 
was ~2000. The lower frame shows the results for the surface 
brightness method together with errors calculated for individ- 
ual pixels. The parameters a and b were either average values 
of all model clouds, including both models with < Ay >= 1.6™ 
and < Ay >- 3.2™ (see Sect. 13. 2> . or were estimated through 
a comparison of surface brightness and extinction maps. The 
colour excess data did not significantly improve the accuracy 
of the surface brightness method. This is not surprising because 
the average parameters were based on simulations with correct 
dust properties and correct value of the background radiation 
field. On the other hand, Figure[2]does show that by using the 
colour excess information similar accuracy can be reached even 
without any such a priori knowledge of the dust properties and 
of the radiation field. Extinction measurements are therefore 
useful to avoid assumptions about dust properties and radiation 
field intensity in the application of our method. The column 
density map obtained with our method still has a much bet- 
ter spatial resolution that the extinction map itself, seven times 
better in this case. 

In the previous discussion, the colour excess information 
was used globally, through the model parameters a and b. 
Alternatively, one could combine the extinction map and the 
column density map obtained from our method, at the lower 
resolution of the extinction map. Spatial variations in the ratio 
of the maps could serve as an indicator of changes in dust prop- 
erties or in the strength of the local radiation field. To test these 
possibilities we modified the model C (< Ay >= 1.6™), first 
re-examining the eff'ect of changing dust properties. The frac- 
tional abundance of dust with Ry = 3.1 was set to depend on 
local density according to the function / = «/(2000 + n). Dust 
with Ry = 5.5 was given relative abundance 1 -/. Column den- 
sities were estimated with both methods (Fig.ll5>. The colour 
excess method gives very accurately the correct average visual 
extinction but, with 2000 simulated background stars, the map 
has significant noise. When the surface brightness method is 
applied using parameters tuned for Ry - 3.1 dust, the column 
densities of the same regions are significantly over-estimated. 
The ratio of the two estimates reflects, therefore, both resolu- 
tion effects and the changes in dust properties. The intensity 
ratios of the NIR bands would, of course, also directly show 
that the assumption of either the dust properties or the radia- 
tion field is not correct. As above, the extinction map is used 
to re-evaluate the average global parameters for the surface 
brightness method. The resulting map is shown in Fig.ll5t. The 
highest column densities are slightly overestimated but the map 
is morphologically extremely accurate. As before, the surface 
brightness data contained observational noise. In Fig.^lall the 
three maps are smoothed to a resolution of six pixels, but the 
signal-to-noise ratios would allow a much higher spatial res- 
olution for the map that was based on the surface brightness 
data. Extinction measurements could also be used in the cre- 
ation of a detailed radiative transfer model of the source (see 
Appendix IaI. That would be particularly useful in the case of 
a strongly anisotropic radiation field. 

Local radiation sources may complicate the derivation of 
the column densities. Depending on the spectrum of these 
sources, it may be difficult to deduce their presence based on 
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Fig. 14. Upper frame: Errors in the simulated extinction map. 
Plots are shown for two versions of model C, with average vi- 
sual extinction 1.6 (circles) or 3.2 magnitudes (squares). The 
resolution of the extinction map corresponds to a FWHM equal 
to seven map pixels. Lower frame: Accuracy of the column 
density estimates based on NIR scattering. Results are shown 
for the two versions of the model C. The parameters are either 
the average values from all models (six clouds, three viewing 
directions and two scalings of average column density; open 
symbols) or were determined by correlating surface bright- 
ness against column densities estimated using the colour excess 
method (solid symbols). 



surface brightness data only. This is true especially if the true 
column density is low and the surface brightness itself is not 
anomalously high. In that case the sources could lead to sig- 
nificant overestimation of the column density. In the following 
we examine how such errors manifest themselves in the es- 
timated column density maps and to what extent the colour 
excess data can be used to correct the results. We added to 
cloud C three point sources, each with a 4000 K black body 
spectrum. The K-band luminosities were 1.5x10''', 3.1x10'^, 
and 6.2x10'^ erg/s/Hz. Figure [T6h-c shows the surface bright- 
ness maps for each source separately, excluding direct radiation 
from the source. The lower frames show again the true column 
densities, the NICER estimates (using 2000 background stars), 
and the estimates based on scattered light, all convolved with a 
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Fig. 15. Model containing variable dust properties, - 3.1 in low density regions and Ry - 5.5 in high density regions. 
The frames are: a. True Ay, b. Ay based on colour excesses of ~2000 background stars, and c. Ay based on scattered surface 
brightness. All maps have been smoothed to a resolution of six map pixels. 



beam with FWHM~6 pixels. In observations, the presence of 
the strongest source would be obvious, based on the presence 
of the point source, high surface brightness values in excess of 
the ISRF intensity, and comparison with the NICER map. The 
second source can still be seen by plotting the difference of the 
two column density estimates. In the surface brightness method 
(using default parameters) the estimated column density at the 
position of source 2 is twice the correct value. However, the ef- 
fect of this source is very local. At the distance of 6 pixels the 
error is still ~IQQ% but falls below 50% at the distance of 12 
pixels. The weakest source is outside dense regions and there- 
fore does not produce A clear peak in the surface brightness 
or estimated column density. There is, however, a diffuse area 
(~1% of the whole map) where the column density is overesti- 
mated by ~3Q%. Because the area has low extinction. Ay < 4, 
the noise is relatively high in the NICER map, and it would be 
difficult to detect the effect of this source in actual observations. 



5. Discussion 

Observations of scattered NIR light provide a means to map 
large areas of interstellar clouds with very high spatial resolu- 
tion. Our simulations showed that one can reach average pixel- 
to-pixel accuracy better than 10% as long as the maximum ex- 
tinction is below Ay ~ 15™. High absolute accuracy requires 
that the parameters of Eq. ^ can be determined with similar 
precision. This requires knowledge of both the radiation field 
and the dust properties. Fortunately, the NIR observations will 
automatically give colour excess data for a large number of 
background stars. Comparison of surface brightness data and 
the colour excess data provides a direct way to determine the 
parameters of Eq. ^ Therefore, the column density estimation 
can be carried out with little a priori information. This is partic- 
ularly important regarding the radiation field because its inten- 
sity is usually unknown and cannot be easily determined. On 



the other hand, the NIR dust properties are generally believed 
to be rather constant among sources of similar type. If the con- 
version between surface brightness and column density is done 
with incorrect parameters, the absolute scaling and the rela- 
tive values, as a function of position or column density, will be 
incorrect. However, the derived column density maps are still 
usually morphologically accurate (see, e.g., Fig.ll5>. 

The use of Eq. [^contains the assumption that the intensity 
of the radiation field and the dust properties are constant within 
the cloud. Spatial variations of dust properties were examined 
in Sect. 14.31 In particular, Fig.ll3lshowed that although varia- 
tion between dust with Ry = 3.1 and Ry =5.5 can introduce 
visible changes in the estimated column density maps, the vari- 
ations are still only at the 10% level. Our simulations showed 
that as long as the clouds are not very optically thick. Ay ~ 15™ 
or below, the errors caused by large scale radiation field gradi- 
ents are ~10% or below. The effect was visible mainly towards 
the edges of relatively homogeneous clouds (see Fig.|8}. 

At large scales the spatial variations in the dust properties 
and the radiation field can be recognized with the help of back- 
ground stars. An extinction map that is based on background 
stars will have a much lower resolution than the surface bright- 
ness maps. Furthermore, because the parameters of Eq. [Qcan 
only be determined using correlations over large areas, rapid 
spatial variations in dust properties may not be taken fully into 
account. Nevertheless, the comparison of extinction and sur- 
face brightness maps will give indication of those variations. 
Some effects can be recognized and corrected based on surface 
brightness data itself. This is discussed further in AppendixlAl 
where we examine the possibility of improving the accuracy of 
the column density estimates through detailed radiative transfer 
modelhng. The effects resulting from radiation field and dust 
property variations are examined further in AppendixiBl 

We have assumed that the observed surface brightness is 
caused entirely by dust scattering. The contribution from emis- 
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Fig. 16. Frames a-c: H-band surface brightness maps for each of the three point sources that were added to model cloud C. Frame 
d-f: Column density estimates for model C containing three point sources. The maps show the true column density, estimates 
based on NICER method and estimates based on scattered surface brightness. Note that the colour scale is logarithmic in the 
upper frames and linear in the lower frames. 



sion lines is small but the role of NIR dust emission remains 
uncertain. At low Ay sightlines a significant part of the surface 
brightness could be caused by emission from very small grains 
or possibly PAHs (Flagey et al. ,2006j . If, as it seems likely, the 
emission is restricted to regions of low extinction and mainly 
lo A 2fim, it will not cause serious problems for the col- 
umn density estimation. If Eq. ^is modified to accommodate 
the additional emission component, colour excess data of back- 
ground stars can still be used to derive the relation between K- 
band surface brightness and column density for higher Ay. This 
question is discussed further in the Appendix lB.4l Further mod- 
ifications to Eq. n will be needed if the surface brightness of 
the background sky is not negligible. At high galactic latitudes 
the extragalactic infrared background is the main background 
component. If completely ignored, it could cause systematic 
errors again at the level of 10% percent (see ApDendix lB.4.3> . 
However, if the parameters of Eq. ^ are determined through 
comparison with extinction data, most of this error is already 
automatically eliminated. 

The scattered light traces the distribution of dust, and the 
total column density can be estimated by assuming a constant 
gas-to-dust ratio. Conversely, the scattered surface brightness 
provides a way to study the relative distribution of gas and 
dust. In the NIR one can easily reach sub-arcsecond resolution. 



which, at the distance of the closest star forming regions cor- 
responds to only ~100AU. Comparison with interferometric 
observations of suitable molecular tracers might reveal entirely 
new phenomena, such as the clustering of dust grains within 
turbulent flows (Padoan et al. l2006> . 

6. Conclusions 

We have examined the use of scattered NIR surface brightness 
as a high resolution tracer of the structure of interstellar clouds. 
Column densities are estimated using the J, H, and K bands. In 
the optically thin case the observed surface brightness is di- 
rectly proportional to column density. In clouds of moderate 
optical thickness. Ay below ~ 15"^°, the saturation of the NIR 
surface brightness can be corrected, and reliable estimation of 
column densities is still possible. We have studied how the col- 
umn density estimation is affected by cloud structure, optical 
depth, and dust property variations. In the numerical models 
the absolute accuracy of the estimated column densities is bet- 
ter than -20%. Large variations in the radiation field or dust 
properties can cause large errors, if not corrected for However, 
even in those cases the errors are mostly systematic, and the 
relative errors between map pixels are significantly lower than 
the systematic errors. 
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Sensitive observations of the surface brightness will auto- 
matically provide colour excess data on background stars that 
can be used as a largely independent measure of the column 
density. We showed that by combining the two methods most 
of the errors caused by anomalous radiation field or dust prop- 
erties can be easily corrected, resulting in a reliable column 
density map with significantly higher resolution than is pos- 
sible with colour excess data alone. Conversely, the compar- 
ison of the two column density estimates can give important 
information on the radiation field and dust property variations 
within the cloud, and even give hints on its three-dimensional 
structure. 
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Appendix A: Use of iterative radiative transfer 
modelling 

The assumption of a uniform radiation field breaks down in optically 
thick regions (see Fig.0)), and even in a homogeneous cloud the in- 
tensity gradients increase the errors of the column density estimates. 
These variations are related to the three-dimensional cloud structure 
and, therefore, can never be fully eliminated. A rough correction is 
possible, if the cloud structure along the line-of-sight can be ignored, 
and only the projected column density map is used. Variations in the 
radiation field can be estimated with radiative transfer calculations. 
An initial model cloud is created according to preliminary estimates of 
the column densities. In the line-of-sight the structure is not restricted 
by observations, and we start by assuming a Gaussian density profile 
with FWHM equal to one third of the cloud size. Radiative transfer 
calculations are used to produce maps of surface brightness, and the 
column densities of the model cloud are corrected iteratively, using the 
ratio of the observed and modelled intensities. The correction is done 
separately for each sightline corresponding to a pixel in the maps. 

The correction procedure was first applied to model F with Ay = 
1.6 and on model D with Ay = 3.2. Cloud F is very inhomogeneous 
and, therefore, a correction without any knowledge of the line-of-sight 
structure is expected to be more difficult than for the more homoge- 
neous cloud D (see Fig.j^. Since the iterative model fitting is some- 
what time-consuming, the model clouds were downsized to 64' cells. 
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Fig.A.l. Comparison of the accuracy of column density es- 
timates obtained with the analytical method (open symbols) 
and with detailed radiative transfer modelling (solid symbols). 
Results are shown for model D with < Ay >= 3.2™ and for 
model F with < Ay >= 1.6™. The plot shows the bias and the 
scatter with respect to the true column densities. 



The relative rms noise of the initial simulated maps was below 1%. 
This applies both to the initial 'observed' surface brightness maps and 
the simulated maps produced during the iterations. 

In Fig. lA. II we compare the resulting column density estimates 
with those obtained using Eq. Q together with parameters a and b 
that were averages of the values obtained from the six models. As 
expected, the modelling works best for the cloud D. In that case esti- 
mates are free from bias and scatter is smaller than obtained with the 
default method. In the latter method, if we used exactly correct param- 
eter values in Eq.Q the bias would be removed but the scatter would 
remain unchanged. For model F the accuracy is rather similar for the 
two methods. The iterative modelling shows a systematic error 2-3%, 
probably as a consequence of the assumption of a smooth line-of-sight 
density distribution. 

The use of Eq. Q includes the effects of a clumpy cloud struc- 
ture since the fit parameters are obtained from inhomogeneous mod- 
els. On the other hand, it ignores any spatial variation on the plane of 
the sky. The iterative modelling worked differently, fitting the spatial 
variations but using a simplistic description for the line-of-sight den- 
sity structure. Remaining errors in the column densities are a direct 
consequence of the true cloud structure differing from the smooth dis- 
tribution used in the modelling. This is more clear in the case of model 
C for which error maps are shown in Fig. lA.2l In this case the iterative 
modelling results in a much larger scatter (see Fig. IA.3> . The maps 
contain intriguing features, where the column density is either sys- 
tematically overestimated or underestimated. The same regions can be 
easily recognized without knowledge of the true column density be- 
cause, when column density is underestimated, the predicted J-band 
intensity from the iterative modelling is relatively higher and the K- 
band intensity lower than in the actual model C. This is caused by 
the fact that, because of the smaller saturation, the K-band intensity 
is more sensitive to changes in the column density. Comparison with 
Fig.|2j) shows that the deviations are not necessarily associated with 
the largest column densities. There could be several explanations for 
such deviations. First, a dense region close to the cloud surface could 
produce more surface brightness (especially in the J-band) and could 



result in an overestimation of the column density. Examination of the 
three-dimensional density distribution does not support this interpre- 
tation. The only region with clearly overestimated column density val- 
ues is located close to the centre of cloud C. Similarly, the strongest 
negative features are caused by filaments that are relatively close to the 
cloud surface. The shape of the condensations may also play a role, 
the surface brightness being slightly higher when the material is dis- 
tributed along a longer sightline. This agrees better with the structures 
seen in the cloud C. When the line-of-sight FWHM of the density dis- 
tribution was used as a further free parameter, the previous fit results 
did indeed improve significantly (see Fig. lA.3b ). 

The scattered surface brightness might also carry some informa- 
tion about the relative positions of the dense regions. If scattering oc- 
curs mostly in the forward direction, the surface brightness samples 
the average radiation field behind the clumps. In that case the clumps 
closest to the observer would have lower ratios between the surface 
brightness and the column density. However, this seems to be a very 
small effect. When model C is observed from the opposite direction, 
the resulting error maps do, within noise, exactly correspond to those 
shown in Fig. lA.2l Finally, the deviations might be caused by the gen- 
eral lack of homogeneity, which would affect the way radiation can 
penetrate the cloud. This seems to be, at least partially, an explanation 
for the underestimated column densities in the cloud centre. This was 
tested in a qualitative way by using the density distribution of cloud D 
in the fitting of model C surface brightnesses. The initial line-of-sight 
densities were random (models D and C are completely independent 
runs), and the column densities corresponding to each map pixel were 
again optimized. The result is shown in Fig. lA.3b . The scatter is now 
rather similar to what was obtained by using Eq.Q The use of an in- 
homogeneous medium eliminates overestimated column densities in 
the central region. However, the negative features seen in Fig.lA^do 
remain. These might be caused by the compactness of the emitting 
regions, as discussed above. Because the deviating regions are visible 
as excess in the simulated J-band maps, it should be possible to make 
a correction for those, for example by using the line-of-sight size of 
the clumps as a further free parameter. Figure lX^ shows results after 
a simple correction where the column density estimates are increased 
by the ratio of modelled and true J-band intensities. This is just to 
illustrate that further improvements are possible. 

When optical depth becomes significant, an increase in the col- 
umn density has gradually less and less effect on the observed inten- 
sity. The saturation depends on the wavelength, and changing inten- 
sity ratios between different bands carry information on the column 
density. In the modelling this can be used to break the degeneracy 
between radiation field and column density. This was tested with the 
model C and repeating the optimization with the scaling of the exter- 
nal radiation field as free parameters. When observational errors were 
included, the correct intensity level could still be recovered with ac- 
curacy better than 10%. This shows that even if the intensity of the 
radiation field were not a priori known, the final uncertainty would 
not be significantly larger than shown in Fig. lA.3l 



Appendix B: Validity of the basic assumptions 

In this section we discuss possible violations of the basic assumptions 
behind the use of Eq.Qand examine their effect on the accuracy of the 
column density estimates. 



B.1. Anisotropy of the radiation field 

The assumption of a uniform radiation field is violated if (1) the back- 
ground radiation is anisotropic, (2) the cloud is sufficiently optically 
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Fig. A.2. Relative errors in the column densities obtained for 
model C, < Ay >- 1.6™, with iterative radiative transfer mod- 
elling. The model assumed Gaussian density distribution along 
the line-of-sight. 



thick to cause significant attenuation in its interior, or (3) the cloud 
contains strong internal radiation sources. The anisotropy of the radia- 
tion field was not studied explicitly. However, it is clear that it becomes 
important only if both conditions (1) and (2) are true. A gradient in the 
intensity of the incoming radiation causes a corresponding gradient in 
the column density estimates. In the models the attenuation of the ra- 
diation field caused errors only at 10% level between cloud edge and 
cloud centre. The sky brightness depends on the Galactic coordinates, 
so that the ISRF is never fully isotropic. The clouds are still illumi- 
nated from all sides, and the Galactic plane fills a large solid angle. 
Therefore, also the difference between cloud edges should exist only 
at a level below ~10%. Although shadowing by dense regions can 
cause variations even at relatively small scales, this can still be iden- 
tified with the help of an extinction map. Once the gradient has been 
identified or the source of the anisotropy is known, it can be corrected, 
either by modelling the radiation field or by varying the fit parameters 
of Eq. Qas the function of position. Direct re-scaling of the derived 
column density estimates does not work as well, because of the non- 
linearity of Eq.Q In the near-infrared the scattering takes place mostly 
in the forward direction. This means that the observed surface bright- 
ness depends more on the sky brightness behind the cloud than in 
other directions, and the scaling between surface brightness and col- 
umn density changes as a function of Galactic coordinates. Of course, 
this can be corrected by taking into account the true sky brightness 
distribution available from DIRBE observations. In the first approx- 
imation the effect is eliminated if the method is re-calibrated using 
background stars. There is a small difference between optically thin 
and optically thick regions, because, in the latter, multiple scattering 
tends to make the field more isotropic. However, when the surface 
brightness is used for the column density estimation, radiation can be 
only moderately optically thick, and the same correction factor should 
be valid for the whole map. 

The radiation field is always attenuated within the cloud. In a ho- 
mogeneous cloud the gradient is systematic and can be estimated by 
modelling (e.g., Fig. lA.H . In practice the radiation field depends on 
the unknown 3D cloud structure. Although an exact correction is not 
possible, our model calculations showed that, by taking into account 
some general properties of the density field, the accuracy of the col- 
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Fig.A.3. Accuracy of column density estimates in the case 
of model C, < Ay >= 1.6™, when three different methods are 
used. Upper frame: Original method using Eq.n(open circles), 
iterative modelling with Gaussian line-of-sight density distri- 
bution (solid circles), and iterative modelling with Gaussian 
line-of-sight density distribution with width taken as a free 
parameter for each sightline (solid triangles). Lower frame: 
Results from iterative modelling with clumpy density distribu- 
tion (squares), and the same after a simple correction based on 
the error in the predicted J-band intensity (triangles). 



umn density estimates can be improved (e.g.. Fig. IA.3> . Our results 
apply to clouds with Ay below 20'". This limit corresponds to a K- 
band optical depth of ~ 1.8. For regions with tk » 1 the scattered 
NIR light can not be used as a tracer of the total column density. As 
the optical depth exceeds unity, the surface brightness saturates and, 
if the region is sufficiently large, finally decreases towards the column 
density maximum. The resulting limb-brightening serves only as an 
indicator of the radius at which the medium becomes optically thick 
for the observed radiation. 



6.2. Internal sources 

Internal radiation sources can cause local increase in the surface 
brightness and decrease the accuracy of the column density estimates 
in active clouds. If a source is deeply embedded its effect on the sur- 
face brightness will be limited to a very small area or it can be readily 
identified, either directly from surface brightness maps or from a com- 
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parison with an extinction map. If the source is outside dense regions, 
its effect may be more subtle, as was seen in Fig. [Tell. A source can 
illuminate distant filaments while, in projection, nearby regions are 
apparently unaffected. As an example, we consider the Chamaeleon I 
cloud where some one hundred embedded sources have been identi- 
fied (Gomez & Kenyon 2001 1. The cloud is located at the distance of 
160 pc, and on the sky the sources are at mean distance of a few arc 
minutes from each other. The brightest source has otk = 10.2"', five 
sources are brighter than 12" but most sources are fainter than 13'". 
Even the strongest source is fainter than the weakest source in Fig.[T6l 
If this source were located in a dense region, it could be visible in the 
surface brightness of scattered light, but the effect would be limited to 
a very small area. Similarly, the contribution of 13'" sources located 
at one arc minute intervals remains almost two orders of magnitude 
below the ISRF intensity. If the source is really embedded inside the 
cloud, the associated elevated NIR dust emission should be limited to 
an even smaller area. Therefore, such embedded sources do not seri- 
ously affect the accuracy of the column density estimates. 

If cloud contains very strong sources, only in very simple cases 
can their effect be fully modelled. However, comparison of extinction 
and surface brightness maps gives a direct way to estimate the effect 
of internal radiation sources and, at the same time, gives clues about 
the 3D structure of the object. When the two column density estimates 
agree, the surface brightness method can be safely used to extend the 
column density estimation to smaller spatial scales. Close to a source, 
the surface brightness maps does contain important information about 
the cloud structure, although weighted according to varying radiation 
field. In practice, the discrepancies between the two column density 
estimates can be caused also by problems in the extinction map, for 
example, if the colour excess data are contaminated by foreground 
stars or background star cluster or galaxies. 

B.3. Dust properties 

Dense regions may be problematic, not only because of the attenuation 
of the radiation field but also because of changes in dust properties. In 
Sect. l4.2.^ we examined how the results are affected if dust properties 
do not agree with the model parameters of Eq. Q Below Ay ~ 10'" 
the effect was only ~20%. This is almost entirely a systematic effect 
and does not affect the relative accuracy between regions of different 
column density. Only above Ay ~ 10'" does the non-linearity of Eq.Q 
produce significant Ay -dependence. In this case no attempt was made 
to correct the values of the erroneous model parameters. In Sect. 14.31 
we examined a more realistic case where dust properties changed ac- 
cording to the local density, from ^y = 3.1 in low density regions to 
Ry = 5.5 at the highest densities. This time a comparison with low- 
resolution extinction map was used to re-estimate the parameters of 
Eq. Q The result was a well calibrated column density map that had 
high spatial resolution and was was morphologically very accurate. 
One set of global parameters was used in Eq. Q Although this did 
work reasonably well for both the dense and the diffuse regions, the 
use of spatially varying parameters might still bring some improve- 
ment. 

Dust properties affect the ratio between scattered surface bright- 
ness and extinction (see Sect. 14. 3t . When grains become larger the 
albedo increases and scatterings are oriented more in the forward di- 
rection. This affects the transport of radiation during the first phase 
(see Sect. |5| and Fig.0 and, for a given extinction, makes the field 
more uniform within the cloud. For clouds with Ay ~ 10'" the vari- 
ations on the plane of the sky are small, and changes in the albedo 
and scattering function can only have a minor effect. The main effect 
comes from the fact that the observed surface brightness is directly 



scaled by the albedo at the position of the last scattering. During the 
second phase the radiative transfer equation is similar as for the ex- 
tinction of the background stars and depends on the extinction cross 
section. The net effect is that the albedo acts as a multiplicative factor 
that scales the surface brightness relative to the extinction. Conversely, 
the ratio of surface brightness and extinction could be used to trace 
variations in the average grain size. However, as shown by Fig. ll5l the 
effect may be small and its detection and separation from all the other 
possible effects is not straightforward. 

B.4. Other sources contributing to the NIR surface 
brightness 

We have assumed that the observed surface brightness is due to dust 
scattering inside the cloud. If there are other significant sources, their 
contribution must be determined and removed or the column den- 
sity estimates are correspondingly in error. We discuss three potential 
factors affecting the NIR surface brightness: Thermal dust emission, 
emission lines, and the role of the NIR background. 

B.4.1 . Thermal dust emission 

At optical wavelengths scattering dominates over the thermal emis- 
sion from cold dust particles. In the mid-infrared the situation is 
already reversed and emission from transiently heated small parti- 
cles and large molecules is the dominant component. In the mid- 
infrared the strong emission bands are generally attributed to PAH 
molecules. However, there is also a broad continuum that extends to 
the NIR, although significantly diminished. We estimate the thermal 
dust emission based on the model of Li & Draine f200H. According 
to their Fig. 8, the expected dust emission in the K band is AI^ ~ 
2 X 10"^' erg/s/sr/H. Excluding all absorption, the expected intensity 
for dust column with Ay = 1.6 is ~ 4.5 x 10"^" erg/s/cm^/sr/Hz. For 
similar Ay the surface brightness of scattered radiation was in our cal- 
culations ~ 2 X 10"" erg/s/cm^/sr/Hz, the number including some re- 
duction caused by absorption. Therefore, in an optically thin case the 
emission could contribute some 20% to the total surface brightness. 
However, while the scattering involves here only one NIR wavelength, 
dust emission depends on the heating that takes place through shorter 
wavelengths. This is especially true for the NIR emission that requires 
grain temperatures of the order of 1000 K and, consequently, heating 
by energetic photons. At Ay = 1'" the attenuation is in the U-band a 
factor of ~5 and even the V-band close to a factor of three. Therefore, 
dust emission is restricted to cloud surfaces, and, even in translucent 
clouds, its contribution to the total surface brightness is only a few per 
cents. In more optically thick clouds the ratio between scattered and 
emitted intensity increases, until also the scattered K-band intensity 
saturates after Ay ~ 10'". In the H-band the expected dust emission 
is only half of the value in the K-band. Moreover, there is evidence 
that the abundance of small grains decreases rapidly in dense clouds 
(e.g., Stepnik et al. .200 3 1, further reducing the relative amount of NIR 
emission expected from optically thick clouds. 

Compared with the continuum, the contribution of possible PAH 
lines and other features is small (e.g., Gordon et al. 2000 Mattioda 
et al. 2005 3). In the Li & Draine model small graphite grains are re- 
sponsible for most of the NIR continuum. For ionized PAHs the NIR 
absorption cross sections may be significantly higher than assumed 
by that model (see Mattioda et al. l2005l i,|2005b). This may affect the 
interpretation of the source of near-infrared continuum but, because 
the modelled small grain contribution was well above the PAH con- 
tinuum, the effect on total predicted emission would be much smaller. 
Furthermore, within most of the volume of our clouds the radiation 
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field is both softer and weaker than the normal ISRF, thus reducing the 
degree of PAH ionization. On the other hand, Sellgren et al. <1996l and 
references therein) observed in several optical reflection nebulae sig- 
nificant NIR dust emission that was similarly attributed to either very 
small, transiently heated dust grains or PAH-type molecules. In the K 
band the emission exceeded the predicted scattering when the temper- 
ature of the central star was larger than ~ 6000 K. One can make a 
conservative assumption that the ratio of the emitted and scattered ra- 
diation follows the ratio of incoming K-band and V-band intensities. 
For a 6000 K black body this ratio is 1.66, while in the case of the 
ISRF (Mathis et al. 1983 ( the ratio is only 0.42, i.e., lower by a factor 
of four. In reality the reduction can be much larger, because the heat- 
ing is likely to require large flux of UV-photons, because of the strong 
attenuation of the optical-UV field in optically thick clouds, and be- 
cause the observed NIR emission may be related to a difi'erence in dust 
properties. Therefore, the results of Sellgren et al. are not in contradic- 
tion with our assumption that in a cloud illuminated by normal ISRF 
the scattering stands for most of the observed surface brightness. 

Further evidence for NIR dust emission comes from DIRBE ob- 
servations. Bernard et al. 1 1994 1 presented a spectrum for the aver- 
age cold diffuse medium at \b\ = 5 deg. After subtraction of zodiacal 
light and stellar emission, the remaining signal was 0.023 MJy/sr at 
3.5/im and 6.7 MJy/sr at lOOyum, corresponding to Ay ~0.5. Only a 
minor fraction of the 3.5/jm signal can be due to scattered light. If 
one assumes that the 3.5 fim value is caused entirely by emission and 
further assumes a flat spectrum, v/y=constant, the predicted K-band 
emission would be 1.5x10"" erg/s/cm"-/sr/Hz. For this low Ay the 
value is higher than the expected intensity of the scattered light, while 
in the J-band the emission would drop to ~10% of the scattered in- 
tensity. Arendt et al. <1998> correlated different DIRBE bands against 
the 100 yum band at \b\ < 30 deg. The correlation coefficient between 
3.5yum and 100/jm intensities was, 1.8x10"^, i.e., lower by almost a 
factor of two. The values might still be overestimated, if the DIRBE 
values contained additional signal from incompletely removed stars or 
from scattering, or were significantly affected by the presence of the 
3.3/im PAH feature. 

Flagey et al. 2006 have recently studied NIR emission using ISO 
and Spitzer data and modelled emission seen at diffuse sightlines. 
Based on their model the extrapolated value at 2;um is ~0.03 MJy sr"' 
for Nu = 10-' cm"^ (see their Fig. 9). At Ay =0.5™ this would cor- 
respond to 1.9x10"" erg/s/cm"^/sr/Hz, i.e., a slightly higher intensity 
than the previous Bernard et al. ( 1994( value. The NIR continuum 
was modelled with a HOOK gray body which is typical of reflection 
nebulae. A lower colour temperature of the continuum would very 
rapidly reduce the predicted K-band intensity. The model of Flagey et 
al. 2006 applied to diffuse medium in the inner Galaxy where, accord- 
ing to their own estimate, the radiation field has intensity three times 
the local ISRF. The emission should be correspondingly smaller in 
local clouds. A critical question is, at what Ay the UV-field has atten- 
uated so much that the scattering dominates also the K-band surface 
brightness. 

The NIR dust emission remains uncertain because of the lack of 
observations and because of the significant uncertainties in the mod- 
elled NIR emission (see also Zubko et al.'2004V There is observational 
evidence that at higher column densities the NIR surface brightness is 
clearly dominated by scattering (e.g., Lehtinen & Mattila 1996 Foster 
& Goodman I2006> . The DIRBE results are consistent with those ob- 
servations, if the dust emission comes mainly from diffuse material. 
In dense clouds, the emission would be restricted to a narrow surface 
layer and it would be distributed rather evenly over the whole cloud. 
Part of the signal may be removed in on-off observations, and the re- 
mainder would not be strongly correlated with the column density. 
Furthermore, if the parameters of Eq.Qare determined from observa- 



tions, e.g., from comparison with extinction data, the contribution of 
dust emission is already taken into account. The functional form of 
Eq.Qcan accommodate some contribution from emission. If emission 
is strong, an additional offset term is required. If the K-band turns out 
to contain significant emission, its use can be restricted to opaque re- 
gions, where dust emission is low. Indeed, the K-band data is essential 
only in very dense regions, where the shorter wavelength data sufi^er 
from significant saturation. Furthermore, the K-band observations are 
relatively more time-consuming, and at low column densities the S/N- 
ratio tends to be worse than for the other bands. Conversely, deep NIR 
observations will enable mapping of the relative intensity of the dust 
emission and scattering. 

B.4.2. Emission lines 

The NIR bands can contain line emission from smaller molecules and 
ions. In cold clouds the excitation is far too low to produce strong NIR 
lines, e.g., from CO and H2 molecules. Significant line emission can 
occur either in photon dominated regions (PDRs) on the cloud sur- 
face or in shocks triggered, for example, by the outflows from young 
stellar objects (YSOs) within the cloud. The presence of strong PDR 
regions would also contradict the assumption of an isotropic back- 
ground. Foster & Goodman 200Qdiscussed the possible contribution 
of H2 lines based on the models of Black & Dishoeck <1987> . They 
noted that the observed J- and H-band intensities showed no indica- 
tion of line emission. Indeed, in the models of Black & Dishoeck, as- 
suming a cloud illuminated with normal ISRF, the effect of H2 should 
remain insignificant. 

YSOs are commonly associated with molecular outflows that are 
driven by fast jets. Along its path the jets cause shocks, Herbig-Haro 
objects, that are visible in the NIR, especially through many H2 lines. 
The strongest line, the H2 1-0 S(l) line at lAlfim, falls within the 
K-band filter. The J and H bands have similarly some contribution 
from, e.g., the excited [Fell] lines at 1.24 and 1.64yum (Reipurth(1999). 
Therefore, along the jet all NIR bands may be contaminated by line 
emission. However, if Herbig-Haro objects are observed, they are usu- 
ally clearly aligned along the jets, which must be masked out from 
subsequent analysis. If no outflow is detected, for example, if even 
the brightest knots fall below the detection limit, the jet will clearly 
have no effects on the column density estimates. At the end of the 
outflow then bowshocks may also cause enhanced NIR surface bright- 
ness. Bowshocks can be much larger than the Herbig-Haro regions, 
and may be located very far from the driving source and the visible 
part of the jet. The shape of a bowshock depends on the density struc- 
ture of the ambient material that is interacting with outflow. Therefore, 
a purely morphological identification of bowshocks may be difficult. 
Strong bowshocks can be recognized based on anomalous intensity 
ratios, while faint ones may go unnoticed and cause some errors in the 
column density estimates. 

B.4.3. Background radiation 

In on-off measurements one measures the difference between the sig- 
nal from the source and the background. If the intensity of the diffuse 
background, 5bg is large behind the source without significantly af- 
fecting the total amount of radiation entering the cloud, the observed 
signal will decrease by an amount of 5bg(l - e^^), where t is the op- 
tical depth of extinction for a sightline through the cloud. As long as 
this term is smaller than the surface brightness due to scattering, the 
dependence between column density and surface brightness will only 
become more shallow. As an empirical description of this dependence 
Eq.Q should remain approximately valid, provided that the parame- 
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ters of the equation take the background into account. Re-calibration 
using the colour excess data would automatically accomplish this. If 
Sbg is equal to the surface brightness of scattered light, column densi- 
ties could no longer be estimated, because the observed signal would 
be zero. If 5bg is larger and sufficiently constant, the cloud structure 
could, of course, be traced by the absorption. 

We will first consider the cosmological IR background (CIRB), 
•Scire- If the absorption due to the Galaxy is ignored, this is an 
isotropic component. If one considers optically thin clouds or a per- 
fectly spherical cloud, the scattered CIRB photons have no effect on 
the observed surface brightness: The distribution of outcoming pho- 
tons is also isotropic. If the cloud is optically thick and inhomo- 
geneous, the radiation field is no longer isotropic inside the cloud. 
Therefore, the amount of radiation scattered away from one line-of- 
sight is no longer equal to the amount of radiation scattered from other 
directions into that direction. However, this imbalance is only a small 
fraction of the scattered CIRB intensity and this variation in the sur- 
face brightness of the scattered CIRB photons can be safely ignored. 
Because of the absorbed radiation the observed intensity decreases by 
Sbg(l - e^^), the term t being the optical depth for absorption. In the 
optically thin case the scattering can be ignored altogether, because 
each scattered photons escapes the cloud and the surface brightness of 
the cloud is identical with the background sky. The depression in the 
surface brightness of optically thicker clouds has been calculated by 
Mattila ( 1976 table A2). When the optical depth is below 4 and the 
contribution from scattering is ignored, the relative error in our esti- 
mate of the change in the surface brightness is less than ~15%. This 
is sufficient accuracy for our purposes. 

The tentative K-band CIRB detection of Gorjian, Wright, and 
Chary CTHt was Sbg = 1.64 x IQ-''' erg/s/cm^/Hz/sr. This high value 
corresponds almost to the mean surface brightness of scattered radi- 
ation in our < Ay >= 1.6™ models. In on-off observations the addi- 
tion of this background would decrease the measured K-band surface 
brightness excess roughly by 5bg(l -e"""*'^"^), the factor 0.049 being 
the ratio of K-band absorption and Ay. The decreased surface bright- 
ness leads to an underestimation of the column densities. The magni- 
tude of this effect is shown in Fig. lB.ll assuming that the parameters of 
Eq.Qare determined without taking the background into account and 
only K-band data are used. If the background level were known or the 
parameters were corrected using the extinction data, the error could be 
almost completely removed. Therefore, in practice it has little effect 
on the accuracy of the column density estimates. The CIRB intensity 
used above is very uncertain, and the true value may turn out to be 
significantly lower (e.g., Dwek, Arendt, Krennrich .2005ii Aharonian 
et al. I2006> . On the other hand, the effect of the CIRB is relatively 
higher in the J- and H-bands, because of the higher optical depth and 
possibly higher CIRB intensity (Hauser & Dwek .20OTl . 

In the modelling we have assumed that the cloud is illuminated 
by the normal ISRF. The diffuse background was not included in 
the calculated surface brightness nor was a background subtraction 
included in the subsequent analysis. If the ISRF is caused by dis- 
crete point sources, no background subtraction is necessary, because 
the surface brightness can be determined on the empty sky between 
sources. However, if the ISRF contains a significant diffuse compo- 
nent, the observed values will again be reduced by 5j||/^""^(l - e"^), 
T being the optical depth for absorption, and qualitatively the effect is 
the same as in Fig. lB.ll If the background intensity is approximately 
isotropic, the main effect comes from absorption, and the scattering 
can again be ignored. Diffuse background could consist of unresolved 
stars or of photons that are scattered from background clouds. In the 
latter case the mapping of the foreground cloud may be impossible, 
because its contribution to the observed surface brightness cannot be 
separated from the background. The effect of faint stars depends on 
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Fig. B.l. Ratio of column density estimates and true column 
density when the analysis ignores an existing background sur- 
face brightness. In the plot the column densities are derived 
based on K-band data and using the average parameters ob- 
tained from the numerical models. The background intensity is 
1.64 X 10"'^ erg/s/cm^/Hz/sr (see text). 

the direction on the sky and the resolution and sensitivity of the ob- 
servations. A crude estimate can be made assuming a slope of 0.35 
in the {mag, ZogioA')-relation (see Sect. I4.3> . The relation is shallow 
enough that, in the total intensity, the contribution of a magnitude bin 
decreases with increasing magnitude. For observations discussed in 
Sect. l4.2l the limiting K-band magnitude is far fainter than 20'". If one 
assumes that the contribution of stars brighter than 20™ can be sep- 
arated from the diffuse background and stars are integrated down to 
40'", the contribution of the unresolved stars remains below 10% of 
the total surface brightness. The net effect on the column density esti- 
mates is again some fraction of this. For a given direction on the sky, 
more accurate estimates and possible correction can be made using 
Galactic models of stellar distribution (e.g., Cohen . 1994,1 . However, 
our estimate is rather conservative, and, in most cases, the effect of 
unresolved stars can be ignored. The Galactic plane may be an im- 
portant exception, because of the anisotropy of the radiation field, and 
the source confusion which makes the separation of the diffuse back- 
ground difficult. 

In the NIR the zodiacal light is by far the strongest diffuse com- 
ponent but, because it has a smooth spatial distribution and it resides 
between the observer and the studied cloud, it has little effect on the 
column density estimation. Some problems can still be caused by sur- 
face brightness gradients across large maps close to the ecliptic plane, 
or, in rare cases, comet trails that introduce small scale structure in the 
surface brightness (e.g., Kelsall et al. ri99Sl Nakamura et al. l2000> . 



